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FOREWORD 


Investigations  of  matrix  operations  have  been  made  in  this  laboratory  for  many 
years.  Subroutines  have  been  prepared  on  the  basis  of  many  methods  of  computation. 
Before  experience  with  the  subroutines  could  be  reported,  the  subroutines  became 
obsolete  as  the  result  of  replacements  of  the  computers.  It  is  the  purpose  of  this 
report  to  document  a  few  surviving  subroutines.  The  manuscript  was  completed  by 
27  October  1977. 


Head,  Strategic  Systems  Department 
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ABSTRACT 


Analysis  and  documentation  are  given  (or  subroutines  which  do  matrix  arithmetic 
and  characteristics  computation.  The  subroutines  make  it  possible  for  the  elements 
of  the  matrix  to  be  in  the  natural  arrangement.  During  inversion  or  trianguloidization 
the  pivot  selection  is  independent  of  the  scaling  of  the  rows  and  columns.  A  subroutine 
does  arithmetic  on  partitioned  matrices. 
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INTRODUCTION 

In  an  early  application,  the  Danilevski  method  was  used  to  compute  the  characteristic 
roots  of  a  large  economic  matrix  on  the  Mark  II  Aiken  Relay  Calculator  for  Professor 
Morgenstern  of  Princeton  University.  Experiments  on  matrix  operations  were  continued 
on  The  Naval  Ordnance  Research  Calculator.  An  early  library  of  subroutines  was  based 
on  a  variety  of  algorithms  which  had  been  extolled  in  the  literature.  Matrix  routines 
which  were  prepared  for  NORC  became  useless  with  the  dismantlement  of  NORC.  They 
were  salvaged  by  a  conversion  to  a  hybrid  version  of  FORTRAN  which  would  run  on 
STRETCH.  Even  the  hybrid  version  of  FORTRAN  became  useless  with  the  dismantlement 
of  STRETCH.  Now  the  computer  in  this  laboratory  is  a  CDC  6700  Computer. 

In  the  meantime  there  have  been  giant  projects  by  large  staffs  of  personnel  for  the 
preparation  of  subroutines  to  perform  all  operations  on  many  forms  of  matrix.  There 
are  the  Math  Sciences  Library  of  Control  Data  Corporation*,  the  System/360  Scientific 
Subroutine  Package  of  the  International  Business  Machines  Corporation2,  and  the 
Subroutine  Library  of  the  International  Mathematical  and  Statistical  Libraries,  Inc.3 
Subroutines  for  the  computation  of  roots  and  vectors  are  available  in  a  Package  of 
Matrix  Eigensystem  Routines*.  Nevertheless,  some  of  the  ideas  in  the  original  NORC 
library  still  are  worth  recognition  and  are  documented  in  the  present  report. 

Available  methods  for  the  transference,  transposition,  addition,  subtraction,  or 
multiplication  of  matrices  are  standard5.  Available  methods  for  the  inversion  of  matrices 
or  for  the  determination  of  characteristic  roots  and  vectors  are  numerous.  The  merits 
of  many  methods  have  been  discussed  by  Wilkinson  *0'12.  The  methods  have  fluctuated 
in  popularity.  When  a  method  was  strongly  defended  it  was  programmed  for  NORC. 
Experience  on  the  computer  has  guided  the  selection  of  methods  for  the  present 
investigation.  Under  consideration  are  only  those  methods  which  are  suitable  for  the 
general  matrix  with  arbitrary  symmetry.  Not  included  are  the  special  methods  which 
are  useful  only  for  sparse  or  symmetric  matrices. 

Matrix  inversion  by  the  evaluation  of  determinants  is  inefficient  for  large  matrices 
because  of  redundancy.  Reduction  of  a  matrix  to  the  coefficients  of  a  characteristic 
polynomial  fails  because  a  polynomial  of  high  degree  cannot  be  evaluated  near  its 
largest  roots. 

The  application  of  sequential  methods  to  asymmetric  matrices  can  get  into  trouble 
from  vanishing  pivots.  Matrix  inversion  by  the  reduction  of  a  matrix  to  diagonal  form 
through  progressive  partitioning  fails  if  the  matrix  has  a  zero  principal  minor 
determinant.  Determination  of  characteristic  roots  by  the  reduction  of  a  matrix  to 
tridiagonal  form  through  progressive  partitioning  fails  if  the  sum  of  products  of 
elements  in  a  row  and  a  column  is  zero. 

Biorthogonalization  methods  have  a  strong  appeal  from  a  philosophical  standpoint 
because  they  work  with  vector  invariants.  Unfortunately,  whenever  two  vector  invariants 
of  different  magnitudes  are  combined  into  a  single  vector,  the  accuracy  of  the  smaller 
vector  is  lost  to  round  off.  Also,  the  biorthogonalization  in  a  three-term  recurrence 
deteriorates  as  the  recurrence  advances,  and  the  deterioration  must  be  kept  in  check 
with  a  multi-term  recurrence. 

Iterative  methods  are  the  most  accurate  insofar  as  they  refer  in  each  cycle  to  the 
original  matrix.  In  too  many  cases  the  rate  of  convergence  of  the  iteration  is  too  slow, 
and  the  convergence  must  be  accelerated  by  various  techniques.  Also,  it  is  necessary 
to  make  a  suitable  choice  of  initial  conditions  for  the  iteration. 
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The  popularity  of  elemental  methods  has  run  full  circle.  In  the  beginning  the  favorite 
method  was  Gauss  elimination  with  pivot  selection.  Now  the  favorite  method  again  has 
become  the  Gauss  elimination  with  pivot  selection. 

The  best  method  for  the  determination  of  the  characteristic  roots  and  vectors 
reduces  the  matrix  to  a  tractable  form  with  the  least  amount  of  damage.  A  symmetric 
matrix  can  be  reduced  to  tridiagonal  form  by  the  Householder  reduction.  An  asymmetric 
matrix  can  be  reduced  to  trianguloid  form  by  permutations  and  eliminations.  The 
roots  of  a  triangular  matrix  are  just  the  diagonal  elements.  Similarity  transformations 
are  applied  in  the  QR  algorithm  to  bring  the  matrix  to  blockwise  triangular  form.  In 
the  present  investigation  the  roots  of  the  trianguloid  matrix  are  derived  by 
Newton -Raphson  iteration. 

Although  there  exists  already  a  subroutine  MATRIX  in  the  CDC  6700  system,  a  new 
subroutine  MTRX  has  been  prepared.  The  new  subroutine  allows  the  interval  between 
columns  to  be  other  than  +1,  and  permits  the  matrices  to  have  their  natural 
arrangement  instead  of  the  transposed  arrangement.  The  subroutine  MATRIX  is  written 
in  machine  language  whereas  the  subroutine  MTRX  is  written  in  FORTRAN.  Timing  trials 
have  shown  that  the  speed  of  MTRX  could  be  increased  if  the  subroutine  were  rewritten 
in  machine  language.  This  would  not  be  true  if  the  FORTRAN  compiler  were  optimum. 

The  NORC  was  a  three-address  binary  coded  decimal  computer.  The  normal  function 
of  a  three-address  instruction  was  arithmetic  on  numbers  in  core  storage.  If  the 
three-address  instruction  were  preceded  by  a  call  to  a  matrix  subroutine,  then  the 
three  addresses  were  reinterpreted  as  specifications  for  whole  matrices.  The  coding 
for  matrix  arithmetic  was  little  more  difficult  than  the  coding  for  number  arithmetic. 

When  matrices  were  too  large  to  fit  into  the  core  memory  of  NORC,  they  were 
partitioned  and  were  stored  on  tape.  The  matrices  could  be  partitioned  rowwise  or 
blockwise  and  they  could  be  prepositioned  or  repositioned  before  operation  The 
computing  loops  were  designed  for  minimum  time  of  operation  and  the  tape  movements 
were  designed  for  minimum  time  of  repositioning.  It  took  six  months  to  design  the 
partitioned  matrix  routine  for  NORC.  Translation  of  the  original  routine  to  FORTRAN 
would  not  be  feasible  because  many  capabilities  which  were  under  program  control 
on  NORC  are  no  longer  under  program  control  in  FORTRAN. 

On  NORC  it  was  possible  to  process  addresses  in  core  with  meta  arithmetic.  In 
FORTRAN  there  is  no  program  control  over  addresses  except  through  indices. 

On  NORC  a  two-dimensional  array  could  be  arranged  in  the  natural  order  of 
mathematics.  In  FORTRAN  the  two-dimensional  array  is  stored  in  transposed  order  in 
memory. 

On  NORC  all  data  were  numeric  and  only  four  of  the  seven  tracks  were  used  on 
tape.  Seven  tape  units  were  available  to  the  programmer.  Each  tape  was  designated 
in  the  program  by  a  symbolic  unit  number.  In  FORTRAN  the  data  are  alphanumeric 
or  binary.  BCD  data  are  written  in  even  parity  while  binary  data  are  written  in  odd 
parity.  Six  tape  units  and  fifty  disk  files  are  available  to  the  programmer.  Disk  files 
are  designated  by  the  same  symbolic  unit  numbers  as  tape  files.  Which  device  actually 
is  associated  with  each  unit  number  is  specified  by  a  device  definition  statement. 

On  NORC  each  block  was  given  a  block  number.  Blocks  could  be  read  forward  or 
backward.  Intermediate  blocks  were  skipped  during  the  search  for  a  specified  block. 
In  FORTRAN  there  is  no  provision  for  record  numbers  and  records  must  be  read  forward 
or  backspaced  one  at  a  time. 

On  NORC  the  writing  of  interrecord  gaps  was  under  program  control.  By  the  use  of 
specially  designed  gaps  between  sentinel  blocks  and  matrix  blocks  it  was  possible  to 
rewrite  a  matrix  anywhere  in  the  interior  of  a  file  of  matrices.  In  FORTRAN  the  record 
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gap  is  not  under  program  control  and  a  matrix  within  a  tape  file  can  be  modified  only 
while  the  entire  file  is  copied  onto  another  tape  file. 

On  NORC  the  writing  of  an  EOF  was  under  program  control.  In  FORTRAN  an  EOF  is 
written  automatically  whenever  a  rewind  or  a  backspace  is  preceded  by  a  write.  A 
rewind  instruction  to  a  tape  unit  causes  delay  while  the  tape  is  rewound  physically, 
whereas  a  rewind  instruction  to  a  disk  unit  merely  resets  an  index. 

For  operations  with  partitioned  matrices  it  is  necessary  to  adopt  limitations  in 
order  to  keep  the  number  of  options  to  a  manageable  level.  Only  conversions  of  format 
are  applied  to  matrices  in  BCD  format.  Transfers  of  data  from  or  to  matrices  in  BCD 
format  are  by  read  and  write  instructions.  Any  number  of  BCD  matrices  may  be  located 
on  each  tape  file.  Matrix  operations  are  performed  only  on  matrices  in  binary  format. 
Transfers  of  data  from  or  to  matrices  in  binary  format  are  by  buffer  in  and  buffer 
out  instructions.  Only  one  binary  matrix  is  assigned  to  each  tape  file. 

On  the  CDC  6700  Computer  the  maximum  order  of  matrix  in  core  under  nominal 
conditions  would  be  144  for  multiplication  and  256  for  inversion.  The  manipulation  of 
larger  matrices  requires  that  the  matrices  be  partitioned  with  the  submatrices  stored 
on  disk.  The  size  of  matrices  is  limited  by  the  availability  of  disk  storage.  The  default 
value  for  the  maximum  number  of  words  in  disk  storage  per  setup  is  896  x  4096  or 
3670016.  The  limit  can  be  increased  fourfold  with  the  aid  of  a  limit  card  to  a  maximum 
of  14680064.  As  long  as  the  maximum  number  of  words  is  not  exceeded,  matrices  may 
be  reread  or  rewritten  indefinitely  on  disk. 


MATRIX  SPECIFICATIONS 

The  elements  of  an  m  x  n  matrix  A  have  the  arrangement 

an-  ■">  ®ln>  "•  aml .  amn 

because  then  they  are  processed  in  serial  order  during  a  matrix-vector  multiplication. 
The  elements  are  stored  in  the  array  AA.  The  matrix  may  have  a  spacing  l  between 
columns  and  a  spacing  fc  between  rows.  The  structure  of  the  matrix  is  specified  in 
core  storage  by  the  words  in  an  array  MA  such  that 

MA(1)=  l  =  interval  between  columns. 

MA(2)  =  to  =  number  of  rows. 

MA(3)  =  n  =  number  of  columns. 

MA(4)  =  k  =  interval  between  rows. 

The  structure  of  the  matrix  is  specified  by  the  four  numbers 

i,  to,  n,  k 

A  compact  vector  of  length  n  may  be  specified  either  as  a  row  matrix  with  the 
specification 

1,  1,  n,  n 

or  as  a  column  matrix  with  the  specification 

1,  n,  1,  1 

The  product  of  a  row  matrix  and  a  column  matrix  is  a  scalar  product  between  two 
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vectors  while  the  product  of  a  column  matrix  and  a  row  matrix  is  the  dyadic  product 
of  two  vectors.  A  compact  matrix  of  size  m  x  n  would  have  the  specification 

1,  m,  n,  n 

while  the  transpose  would  have  the  specification 

n,  n,  m,  1 

for  the  same  actual  arrangement.  A  submatrix  in  a  larger  matrix  of  order  N  would 
have  the  specification 

1,  m,  n,  N 

A  null  matrix  with  the  specification 

1,  n,  n,  n 

can  be  created  by  transfer  from  a  matrix  which  has  the  specification 

0,  n,  n,  0 

and  has  a  single  element  equal  to  zero.  The  null  matrix  is  converted  into  the  identity 
matrix  when  the  diagonal  elements  with  the  specification 

1,  n,  1,  n  +  1 

are  created  by  transfer  from  a  matrix  with  a  single  element  equal  to  unity.  The  real 
parts  and  the  imaginary  parts  of  a  matrix  of  complex  numbers  would  be  matrices 
with  the  specification 

2,  m,  n,  2n 

All  matrices  which  are  involved  in  a  matrix  operation  must  have  compatible  dimensions. 

All  matrices  on  tape  files  are  compact.  The  structure  of  a  matrix  is  specified  on 
tape  files  by  a  sentinel  block  which  precedes  the  matrix  block  The  sentinel  block 
specifies  the  dimensions  m,  n  of  the  matrix.  In  the  case  of  a  BCD  matrix  the  dimensions 
are  punched  on  a  single  card  which  precedes  the  cards  of  the  matrix  In  the  case  of 
a  binary  matrix  the  dimensions  are  written  in  a  2-word  record  which  precedes  the 
77i7i-word  record  of  the  matrix. 


MATRIX  NOTATION 

For  definitions  and  theorems,  reference  may  be  made  to  many  texts  on  matrix 
algebra.  An  m  x  n  matrix  A  is  an  ordered  rectangular  array  of  numbers  a with  the 
following  arrangement 
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Matrix  A  may  be  transposed  by  interchange  of  rows  and  columns  into  matrix  A'.  Every 
nonsingular  square  matrix  A  has  an  inverse  matrix  A-1  such  that 

A  ■  A-1  =  I  (2) 

where  I  is  the  identity  matrix.  The  inverse  matrix  is  given  in  accordance  with  Cramer’s 
rule  by  the  equation 


^  In 

HI  HI  U| 

4l 

HI  HI  ■"  HI 

^«1  ^nj  ^nn 

HI  ”■  HI  ■“  HI 


(3) 


where  the  determinant  AiS  is  the  cofactor  of  the  element  a}i  in  the  determinant  |A|. 
The  Cramer’s  rule  is  the  rule  to  use  for  matrices  up  to  the  third  order,  but  it  requires 
too  much  redundant  computation  for  matrices  of  larger  order. 


MATRIX  ARITHMETIC 

Analysis 

Operations 

When  matrix  A  is  transferred  to  matrix  B  the  elements  are  related  by  the 
transformation 

btj  -*  ai}  ( 1  5  i  i  m,  1  S  ;  §  n)  (4) 

When  matrix  A  is  transposed  into  matrix  B  the  elements  are  related  by  the 
transformation 

bt/  -*  ati  (l  §  i  ^  n,  1  S  j  i  m)  (5) 

If  the  matrix  C  is  the  sum  A  +  B  of  the  matrices  A,  B,  the  elements  are  related  by  the 
transformation 

Cy  -*  aif  +  bif  (1  S  i  g  m,  1  S  j  §  n)  (6) 

and  if  the  matrix  C  is  the  difference  A  -  B  of  the  matrices  A,  B,  the  elements  are 
related  by  the  transformation 

cif  -»  at>  -  6y  (1  £  i  S  m,  1  §  j  g  n)  (7) 

If  the  matrix  C  is  the  product  A  •  B  of  the  matrices  A,  B,  the  elements  are  related  by 
the  transformation 

m 

ctt=E  (IStSf.  1  SfcSn)  (8) 

During  the  inversion  of  a  matrix  A,  it  is  multiplied  by  a  sequence  of  matrices  which 
convert  matrix  A  into  matrix  I  and  convert  matrix  I  into  matrix  A-1.  The  matrix  A  is 
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replaced  by  a  matrix  B  which  evolves  into  the  matrix  A-1.  When  the  element  bT,  is 
used  as  the  pivot  for  elimination  the  elements  are  replaced  in  accordance  with  the 
transformations 


.  .  bia  , 

ba  -  bn  -  brj 


bi, 

b « -  -  f- 


bij  .  +  ~ 


(i  *  r,j  /  s)  (9) 
(i  *  r,j  =  s)  (10) 
(i  =  r,j  *  s)  (11) 
(i  =  r,;  =  s)  (12) 


In  the  elimination  process  the  ith  row  is  replaced  by  the  sum  of  the  ith  row  and  a 
fraction  e  of  the  rth  row  where  e  is  given  by  the  equation 


Let  the  partially  reduced  matrix  be  expressed  by  the  equation 
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e  of  the  rth  row.  When  the  inverse 
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is  used  as  a  postfactor,  the  rth  column  is  replaced  by  the  difference  of  the  rth  column 
and  the  fraction  e  of  the  ith  column.  Once  the  elements  in  all  rows  of  a  column  except 
the  pivotal  row  are  eliminated,  they  remain  undisturbed  while  the  elements  in  all  rows 
of  another  column  with  a  different  pivotal  row  are  eliminated. 

During  elimination  the  indices  r,  s  are  stored  in  arrays.  In  the  index  array  R  the 
indices  r  are  stored  in  ascending  order  of  the  indices  s,  while  in  the  index  array  S 
the  indices  s  are  stored  in  ascending  order  of  the  indices  r. 

When  the  indices  r,  s  have  been  selected  for  the  pivot,  then  the  rth  row  and  the 
sth  column  are  locked  out  of  any  further  consideration  as  a  source  of  pivots. 

If  the  pivots  are  not  diagonal,  the  product  T  of  the  eliminants  transforms  the  matrix 
A  into  a  permutation  P,  and  transforms  the  permutation  P  into  the  matrix  B  in 
accordance  with  the  equations 

P  =  T  A  B  =  T  ■  P  (17) 

The  elimination  of  matrix  T  leads  to  the  equation 

A  1  =  P'  B  P’  (IB) 

where  the  permutation  satisfies  the  equation 

Pl  =  P’  (19) 

The  permutation  is  removed  by  a  succession  of  interchanges  between  rows  and  between 
columns. 

The  values  of  r  for  two  rows  to  be  interchanged  are  derived  from  the  array  R.  The 
rows  are  interchanged  and  the  array  is  adjusted  until  the  values  of  r  are  in  ascending 
order.  The  values  of  s  for  two  columns  to  be  interchanged  are  derived  from  the  array  S. 
The  columns  are  interchanged  and  the  array  is  adjusted  until  the  values  of  s  are  in 
ascending  order. 

Matrix  arithmetic  with  real  matrices  can  be  adapted  to  matrix  arithmetic  with 
complex  matrices.  Let  complex  matrices  be  given  by  the  expressions 

A  +  iB  C  +  iD  (20) 

Then  complex  transfer  is  expressed  by  the  substitutions 

C  -» A  D  ->  B  (21) 

Complex  addition  is  expressed  by  the  equation 

(A  +  iB)  +  (C  +  iD)  =  (A  +  C)  +  (B  +  D)  i  (22) 
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Complex  subtraction  is  expressed  by  the  equation 

(A  +  i  B)  -  (C  +  i  D)  =  (A  -  C)  +  (B  -  D)  i  (23) 

Complex  multiplication  is  expressed  by  the  equation 

(A  +  i  B)  •  (C  +  i  D)  =  (A  •  C  -  B  D)  +  (A  D  +  B  ■  C)  i  (24) 

Complex  inversion  is  expressed  by  the  equation 

(A  +  iB)"‘  =  (A  +  B  •  A'1  •  B)-‘  -  (B  +  A  ■  B~‘  •  A)-1  i  (25) 

Thus  complex  matrix  arithmetic  can  be  achieved  with  the  aid  of  a  subroutine  for  real 
matrix  arithmetic.  However,  complex  inversion  with  real  matrices  is  less  efficient  than 
direct  inversion  of  complex  matrices. 

Pivot  Selection 

In  accordance  with  Cramer's  rule,  each  element  occurs  in  all  cofactors  except  those 
for  the  same  row  or  column.  The  contribution  of  the  element  to  the  inverse  is  the 
product  of  the  element  and  its  cofactor  in  each  determinant  in  which  the  element 
appears.  The  relative  accuracy  with  which  the  element  must  be  preserved  in  the 
inversion  is  relaxed  when  either  the  element  or  its  cofactor  is  especially  small.  An 
optimal  pivot  selection  would  preserve  the  element  with  just  whatever  accuracy  is 
required  for  the  final  inverse.  The  required  accuracy  depends  not  only  on  the  nature 
of  the  matrix  but  also  on  the  use  of  the  inverse.  Without  prior  information  for  each 
case  it  is  necessary  to  adopt  a  general  policy  of  pivot  selection. 

Various  tactics  have  been  used  for  pivot  selection.  Three  options  are  provided  by 
the  subroutine  MATRIX.  In  the  nonpivoting  option  the  pivots  are  the  diagonal  elements 
in  succession.  In  the  partial  pivoting  option  the  pivots  are  the  largest  elements  in 
each  column  in  succession.  In  the  complete  pivoting  option  the  pivots  are  the  largest 
elements  in  the  remainder  of  the  matrix. 

Nonpivoting  is  sufficient  for  a  matrix  with  diagonal  dominance.  Partial  pivoting  is 
inaccurate  for  a  special  matrix  which  has  been  discovered  by  Wilkinson10.  Successive 
eliminations  have  the  effect  of  doubling  the  size  of  elements  in  the  last  row  and 
column.  In  the  last  step  the  subtraction  of  the  large  elements  reduces  them  to  zero 
or  rounding  error.  Permutation  of  the  columns  of  the  matrix  converts  it  into  a 
Hessenberg  matrix  to  which  partial  pivoting  can  be  applied  with  success.  Experience 
with  many  matrices  has  led  Wilkinson  to  favor  complete  pivoting. 

The  relative  rounding  error  in  the  elimination  process  is  independent  of  any  scaling 
of  the  rows  or  columns  which  does  not  change  the  sequence  of  pivots.  If  the  largest 
element  were  taken  as  pivot,  then  the  sequence  of  pivots  could  be  upset  if  the  selected 
element  were  diminished  by  a  selective  scaling.  The  matrix  should  be  prescaled  before 
pivot  selection,  or  the  pivot  selection  should  be  independent  of  scaling. 

If  the  element  6r,  is  the  pivot  for  the  modification  of  the  element  by,  then  the 
accuracy  of  the  element  by  is  retained  to  within  rounding  error  if  the  elements  satisfy 
the  criterion 
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On  the  other  hand,  if  br}  were  the  pivot,  then  the  accuracy  of  the  element  bu  would 
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be  retained  to  within  rounding  error  if  the  elements  satisfied  the  criterion 
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Whichever  quotient  exceeds  unity  by  the  least  amount  would  point  to  the  better  pivot. 
A  symmetric  criterion  would  be  given  by  the  difference  between  the  quotients. 
Subtraction  of  the  two  quotients  gives  a  difference  which  ranges  from  -<*>  to  +■>  and 
gives  too  much  emphasis  on  the  smallest  elements  where  the  requirements  for  relative 
accuracy  can  be  relaxed.  Less  emphasis  on  the  smallest  elements  is  given  by  the 
criterion 
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Then  the  quotients  range  in  value  from  -1  to  +1  and  are  independent  of  scaling. 
Insofar  as  all  of  the  quotients  for  one  pivot  are  superior  to  all  of  the  quotients  for 
the  other  pivot,  it  is  necessary  only  to  determine  the  sign  of  the  criterion  in  order 
to  determine  which  is  the  better  pivot. 


Programming 

Matrix  arithmetic  is  executed  by  reference  to  the  following  subroutine. 

SUBROUTINE  MTRX  (MO.  AA,  MA,  AB,  MB.  AC,  MC) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  MATRIX  ARITHMETIC 

•A******************************************************************************************* 

The  mode  of  operation  is  given  in  argument  MO.  The  matrices  A,  B,  C  are  given  or 
stored  in  the  arrays  AA,  AB,  AC.  The  specifications  of  the  matrices  are  given  in  the 
arrays  MA,  MB,  MC.  The  specification  for  a  matrix  AX  consists  of  the  array  MX  as  given 
in  the  following  table. 


Address 

Specification 

MX(  1 ) 

Interval  between  columns 

MX(2) 

Number  of  rows 

MX(3) 

Number  of  columns 

MX(4) 

Interval  between  rows 

of 

operations  and  call  lines  is  given 

in 

the 

following  table. 

Operation 

Call  : 

Line 

B 

=  A 

CALL 

MTRX 

(0. 

AA, 

MA, 

AB, 

MB) 

B 

=  A' 

CALL 

MTRX 

(1. 

AA, 

MA, 

AB. 

MB) 

C 

=  A  +  B 

CALL 

MTRX 

(2, 

AA. 

MA, 

AB. 

MB,  AC,  MC) 

C 

=  A  -  B 

CALL 

MTRX 

(3, 

AA. 

MA, 

AB. 

MB.  AC.  MC) 

C 

=  A  B 

CALL 

MTRX 

(4, 

AA. 

MA, 

AB, 

MB,  AC.  MC) 

B 

=  A’1 

CALL 

MTRX 

(5. 

AA, 

MA, 

DA) 

During  matrix  inversion  the  matrix  is  replaced  by  its  inverse  and  its  determinant  is 
stored  in  the  argument  DA. 


PARTITIONED  MATRIX  ARITHMETIC 


Analysis 

The  matrices  are  partitioned  into  submatrices  by  straight  lines  parallel  to  rows  or 
columns.  The  numbers  of  rows  and  columns  of  submatrices  are  indicated  by  a  sentinel 
record  which  precedes  the  whole  matrix,  while  the  numbers  of  rows  and  columns  of 
elements  in  each  submatrix  are  indicated  by  a  sentinel  record  which  precedes  the 
submatrix.  All  matrices  which  are  involved  in  a  matrix  operation  must  have  compatible 
partitioning.  Matrices  which  are  to  be  inverted  must  have  square  diagonal  submatrices. 
All  matrices  must  be  in  the  format  which  is  compatible  with  binary  buffer  instructions. 

The  matrix  A  is  an  ordered  rectangular  array  of  submatrices  AXl/.  When  matrix  A  is 
transferred  to  matrix  B  the  submatrices  are  related  by  the  transformation 

BXv  -*  Ax„  (1  SXSm,  1  sygn)  (29) 

When  matrix  A  is  transposed  into  matrix  B  the  submatrices  are  related  by  the 
transformation 

BXv-Aix  (1  §\§n,  1  gi/Sm)  (30) 

If  the  matrix  C  is  the  sum  A  +  B  of  the  matrices  A,  B.  the  submatrices  are  related  by 
the  transformation 

Cx„  - Ax„  +  BXv  (1  SASm,  1  Si/Sn)  (31) 

and  if  the  matrix  C  is  the  difference  A  -  B  of  the  matrices  A,  B,  the  submatrices  are 
related  by  the  transformation 

CXv  AXv  -  BXv  (iasm,lSi/Sn)  (32) 

If  the  matrix  C  is  the  product  A  B  of  the  matrices  A,  B,  the  submatrices  are  related 
by  the  transformation 

m 

CXv-*  Y,  A^  B„V  (1  SASf.  1  Si/Sft)  (33) 

M-l 

When  matrix  A  is  inverted  into  matrix  B  by  elimination  with  diagonal  pivots  the 
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submatrices  are  modified  in  accordance  with  the  transformations 


~  BMV 

(A  #  fj.,  i/  *  n)  (34) 

-  ~  B^ 

(\  *  H,  v  =  n)  (35) 

B%v  ■*  +  B^  Buv 

(A  =  fi.  v  *  n)  (36) 

BXv  ' 

(A  =  n,  v  =  m)  (37) 

Submatrices  are  read  or  written  from  tape  files.  Arithmetic  operations  on  submatrices 

are  executed  in  core  storage. 

Programming 

The  matrices  in  tape  flies  are  read  into  the  core  memory  and  are  written  from  the 
core  memory  by  references  to  subroutines  Whether  the  matrices  are  in  BCD  format 
or  are  in  binary  format  is  determined  by  whether  a  control  index  is  zero  or  unity 
Submatrices  are  processed  sequentially  or  are  rewound  in  order  to  eliminate  wasteful 
backspacing.  The  programming  must  be  planned  so  that  submatrices  become  available 
when  they  are  needed.  Submatrices  which  are  available  but  will  not  be  needed  until 
later  must  be  stored  in  temporary  files.  Whenever  submatrices  are  processed  more 
than  once  they  must  be  stored  on  a  pair  of  tape  files  by  references  to  subroutines 
Which  tape  file  is  read  or  written  is  determined  by  whether  a  control  index  is  odd  or 
even. 

Partitioned  matrix  arithmetic  is  executed  by  reference  to  the  following  subroutine 

SUBROUTINE  PMTX  (MO.  NA,  AA,  MA,  NB.  AB.  MB.  NC.  AC.  MC.  NP.  NO) 

FORTRAN  SUBROUTINE  FOR  PARTITIONED  MATRIX  ARITHMETIC 


The  mode  of  operation  is  given  in  argument  MO.  The  tape  file  numbers  for  matrices 
A.  B,  C  are  given  in  the  arguments  NA,  NB,  NC.  The  submatrices  are  stored  temporarily 
in  the  arrays  AA,  AB,  AC.  The  specifications  for  submatrices  are  stored  in  the  arrays 
MA,  MB,  MC.  The  tape  file  numbers  for  temporary  storage  are  given  in  the  arguments 
NP,  NQ.  Calls  are  made  to  Subroutine  MTRX,  Subroutine  RDMTRX,  Subroutine  WRMrRx, 
Subroutine  RDMXDX,  and  Subroutine  WRMXDX. 

For  the  transference  B  =  A  the  call  line  to  the  subroutine  is 

CAl  1  PMTX  (0,  NA,  AA,  MA.  NB,  AB.  MB) 

The  submatrices  are  read  from  tape  file  NA  into  array  AA,  are  transferred  from  array 
AA  to  array  AB,  and  are  written  from  array  AB  onto  tape  file  NB. 

For  the  transposition  B  -  A'  the  call  line  to  the  subroutine  is 

CALL  PMTX  ( 1 .  NA,  AA,  MA,  NB,  AB.  MB) 

The  whole  matrix  is  read  as  many  times  as  it  has  columns.  The  submatrices  are  read 
from  tape  file  NA  into  array  AA,  the  submatrices  in  a  column  are  selected,  and  the 
transposed  submatrices  are  written  from  array  AB  onto  tape  file  NB 
For  the  addition  C  =  A  +  B  the  call  line  to  the  subroutine  is 

CALL  PMTX  (2,  NA,  AA,  MA,  NB,  AB,  MB,  NC,  AC,  MC) 

The  submatrices  are  read  from  tape  files  NA,  NB  into  arrays  AA,  AB,  and  the  sum  is 
written  from  array  AC  onto  tape  file  NC. 


For  the  subtraction  C  =  A  -  B  the  call  line  to  the  subroutine  is 

CALL  PMTX  (3.  NA.  AA.  MA.  NB.  AB.  MB.  NC.  AC.  MC) 

The  submatrices  are  read  from  tape  flies  NA.  NB  into  arrays  AA,  AB.  and  the  difference 
is  written  from  array  AC  onto  tape  file  NC. 

For  the  multiplication  C  =  A  •  B  the  call  line  to  the  subroutine  is 

CALL  PMTX  (4,  NA,  AA.  MA.  NB,  AB.  MB,  NC,  AC,  MC,  NP,  NO) 

The  prefactor  A  is  read  once  while  the  postfactor  B  is  read  as  many  times  as  there 
are  rows  in  the  prefactor  A.  The  submatrices  are  read  from  tape  flies  NA,  NB  into 
arrays  AA,  AB.  The  product  of  the  submatrices  is  stored  in  array  AC.  If  the  number  of 
columns  of  the  prefactor  is  greater  than  one,  the  accumulation  of  the  submatrices 
for  one  row  are  stored  temporarily  on  tape  files  NP,  NQ.  When  the  accumulation  of  the 
submatrices  is  complete,  the  submatrices  for  one  row  are  written  from  array  AC  onto 
tape  file  NC. 

For  the  inversion  B  =  A-1  the  call  line  to  the  subroutine  is 

CALL  PMTX  (5,  NA,  AA.  MA,  NB,  AB,  MB,  NC,  AC,  MC.  NP,  NQ) 

The  matrix  to  be  inverted  is  located  initially  on  tape  file  NA.  During  transformation 
by  Jordan  elimination  the  successive  versions  of  the  matrix  are  stored  alternately  on 
the  pair  of  tape  files  NA,  NB.  At  the  completion  of  computation  the  inverse  matrix  is 
stored  on  both  tape  files  NA,  NB.  At  the  beginning  of  each  cycle  of  elimination  the 
submatrices  in  the  pivotal  row  are  written  on  whichever  of  the  tape  files  NP,  NQ  is 
unoccupied.  The  pivotal  submatrix  is  selected  and  its  inverse  is  stored  in  array  AA. 
The  submatrices  in  the  pivotal  row  are  read  back  into  array  AB,  are  multiplied  by  the 
inverse  in  array  AA,  and  the  products  are  written  from  array  AC  onto  the  tape  file  NC. 
The  pivotal  submalrix  on  tape  file  NC  is  replaced  by  the  inverse  of  the  pivotal  submatrix. 
Before  the  transformation  of  each  row  of  submatrices  the  submatrix  in  the  pivotal 
column  is  read  initially  from  the  first  column  of  the  original  matrix,  or  is  read 
subsequently  from  one  of  the  tape  files  NP,  NQ  into  array  AA  For  each  elimination  a 
submatrix  is  read  from  tape  file  NC  into  array  AB,  is  multiplied  by  the  submatrix  in 
array  AA,  and  the  product  is  stored  in  array  AC.  The  submatrix  to  be  modified  is  read 
into  array  AB  from  one  of  the  tape  files  NA,  NB,  and  after  modification  it  is  written 
back  out  from  array  AB  onto  the  other  of  the  tape  files  NA,  NB.  In  the  pivotal  column 
the  submatrix  to  be  modified  is  replaced  by  the  difference  between  zero  and  the 
product  in  array  AC.  In  the  next  column  after  the  pivotal  column  the  modified  submatrix 
is  written  also  on  the  other  of  the  tape  files  NP,  NQ.  At  the  end  of  each  cycle  of 
elimination  the  submatrices  are  read  from  tape  file  NC  and  are  written  on  the  other 
of  the  tape  files  NA,  NB. 

A  sentinel  record  and  a  matrix  record  are  read  from  a  tape  file  into  an  array  by 
reference  to  the  following  subroutine 

SUBROUTINE  RDMTRX  (MO,  NA,  AA,  MA) 

***#*•****•*****••**•*****•«•****•*****••*«•*•*•************«**••*♦*«*********♦*•*****«****** 
FORTRAN  SUBROUTINE  TO  READ  MATRIX  RECORD 


The  mode  of  operation  is  given  in  argument  MO.  The  mode  of  operation  is  zero  for 
BCD  format  and  is  unity  for  binary  format.  The  tape  file  number  is  given  in  argument  NA, 
the  matrix  is  stored  in  the  array  AA,  and  the  specification  is  stored  in  the  array  MA. 
A  sentinel  record  and  a  matrix  record  are  written  on  a  tape  file  from  an  array  by 
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reference  to  the  following  subroutine. 

SUBROUTINE  WRMTRX  (MO.  na.  aa,  MA) 


FORTRAN  SUBROUTINE  TO  WRITE  MATRIX  RECORD 
*#****★******♦*******♦♦***♦♦♦*♦**♦♦**#**♦*♦♦******•***••**#*♦*♦***<*****•****♦**##•#*♦**♦;*#** 

The  mode  of  operation  is  given  in  argument  MO.  The  mode  of  operation  is  zero  for 
BCD  format  and  is  unity  for  binary  format.  The  tape  file  number  is  given  in  argument  NA, 
the  matrix  is  given  in  the  array  AA.  and  the  specification  is  given  in  the  array  MA. 

A  sentinel  record  and  a  matrix  record  are  read  from  a  duplex  file  into  an  array 
by  reference  to  the  following  subroutine. 

SUBROUTINE  RDMXDX  (MO,  NA,  NB,  AA,  MA) 


FORTRAN  SUBROUTINE  TO  READ  MATRIX  FROM  DUPLEX  FILE 


The  mode  of  operation  is  given  in  argument  MO.  The  tape  file  numbers  of  a  pair  of 
files  are  given  in  the  arguments  NA,  NB,  the  matrix  is  stored  in  the  array  AA,  and  the 
specification  is  stored  in  the  array  MA.  The  matrix  is  read  from  tape  file  NA  if  the 
mode  of  operation  MO  is  odd,  but  the  matrix  is  read  from  tape  file  NB  if  the  mode  of 
operation  MO  is  even. 

A  sentinel  record  and  a  matrix  record  are  written  on  a  duplex  file  from  an  array 
by  reference  to  the  following  subroutine. 

SUBROUTINE  WRMXDX  (MO,  NA.  NB,  AA,  MA) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  TO  WRITE  MATRIX  ON  DUPLEX  FILE 
#*»#*****•»*••**•****•*«•**«••**«««*«•***!«•*«**••«•***********»«**********•»**•**•*********«* 

The  mode  of  operation  is  given  in  argument  MO  The  tape  file  numbers  of  a  pair  of 
files  are  given  in  the  arguments  NA,  NB,  the  matrix  is  given  in  the  array  AA,  and  the 
specification  is  given  in  the  array  MA.  The  matrix  is  written  on  tape  file  NA  if  the  mode 
of  operation  MO  is  odd.  but  the  matrix  is  written  on  tape  file  NB  if  the  mode  of 
operation  MO  is  even. 


CHARACTERISTIC  ROOTS  AND  VECTORS 

Analysis 
Regular  Matrix 

A  covariant  characteristic  vector  aM  of  a  matrix  A  satisfies  the  equation 

(A  -  cc„  I)  a„  =  0  (38) 

where  is  the  /ith  characteristic  root.  A  contravariant  characteristic  vector  av  of 
the  matrix  A  satisfies  the  equation 

av  •  (A  -  a„  I)  =  0  (39) 

where  av  is  the  i/th  characteristic  root.  That  av  and  aM  are  orthogonal  if  a„  is  not 
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identical  with  follows  from  the  identity 

a"  •  A  •  a„  =  a„a1'  •  a„  =  alt*v  a„  (40) 

If  the  characteristic  roots  all  are  distinct  then  the  characteristic  vectors  are  distinct 
and  can  be  normalized  in  accordance  with  the  equation 

a*'  a„=l  (41) 

The  matrix  of  order  n  can  be  synthesized  from  idempotents  in  accordance  with  the 
equation 

A  =  £  «„a„a1'  (42) 

i/-i 

If  one  or  more  of  the  characteristic  roots  are  zero  the  matrix  is  singular  If  the 
characteristic  roots  are  equal  for  two  or  more  distinct  vectors,  then  any  linear 
combination  of  the  vectors  is  a  characteristic  vector.  Two  or  more  characteristic  roots 
may  be  equal  because  their  characteristic  vectors  are  collinear,  in  which  case  the 
matrix  is  defective. 

Nonzero  characteristic  vectors  are  possible  only  if  the  characteristic  roots  are 
solutions  of  the  determinantal  equation 

p(ot)  =  |A  -  a  I|  =  0  (43) 

A  direct  expansion  of  the  characteristic  determinant  of  nth  order  would  lead  to  a 
characteristic  polynomial  of  the  nth  degree.  A  direct  expansion  is  feasible  for  matrices 
up  to  the  third  order,  but  larger  matrices  must  be  reduced  to  a  simpler  form.  A 
prefactor  S  and  a  postfactor  S“‘  may  be  applied  to  the  matrix  A  to  convert  it  into  a 
matrix  B  in  accordance  with  the  similarity  transformation 

B  =  S  •  A  •  S"‘  (44) 

Then  the  characteristic  vector  aM  of  matrix  A  is  related  to  the  characteristic  vector  bM 
of  matrix  B  by  the  equation 

aM  =  S-  (45) 

and  the  characteristic  vector  aw  of  matrix  A  is  related  to  the  characteristic  vector  b* 
of  matrix  B  by  the  equation 

a1'  =  bv  S  (46) 

The  characteristic  roots  of  matrix  A  are  undisturbed  by  a  similarity  transformation 
into  matrix  B. 

Let  the  components  of  the  covariant  vector  aM  with  respect  to  a  set  of  base  vectors  be 


(47) 


and  let  the  components  of  the  contravariant  vector  a1'  with  respect  to  the  set  of  base 
vectors  be 

a1'  =  K.  <«  (48) 
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Orthogonality  between  the  characteristic  vectors  is  expressed  by  the  equation 

aM  a "  =  £  a*ai  =  6 »  (49) 

*«1 

where  <5*  is  zero  when  n  *  v  but  unity  when  n  =  v.  The  component  a*  is  the  element 
in  the  ith  row  and  the  /xth  column  of  a  modal  matrix  V,  and  the  component  aj'  is  the 
element  in  the  7th  column  of  the  i/th  row  of  the  inverse  matrix  V~l.  The  similarity 
transformation 

V  ‘  A  V  (50) 


reduces  the  matrix  A  to  diagonal  form. 

Trianguloid  Matrix 

A  Hessenberg  or  trianguloid  matrix  has  nonzero  elements  on  or  below  the  diagonal 
next  above  the  main  diagonal,  and  has  zero  elements  elsewhere  above  the  upper 
diagonal. 

The  characteristic  determinant  is  evaluated  for  any  value  of  a  by  a  progressive 
recurrence.  Let  Am  be  the  value  of  the  mth  principal  minor  in  the  upper  left  corner 
of  the  matrix.  The  mth  principal  minor  is  expanded  into  the  sum  of  products  of 
elements  in  the  mth  column  and  their  cofactors.  The  principal  minors  are  generated 
by  a  recurrence  which  starts  with  the  equation 

A,  =  a„  -  a  (51) 

and  continues  with  the  equation 

~  (Omm  ”  ~  (52) 

where  6m  ,  is  that  cofactor  which  is  obtained  by  the  permutation  of  the  (m  -  l)th 
row  and  the  mth  row  of  the  mth  principal  minor.  The  cofactors  are  generated  by  a 
recurrence  which  starts  with  the  equation 

<5,  =  <V>  (53) 

and  continues  with  the  equation 

<5*  "  amk^k  1  ~  ak-\,k&k-l  (54) 

until  k  -  m  1  The  first  derivatives  of  the  principal  minors  are  generated  by  a 
recurrence  which  starts  with  the  equation 


ciA, 

da 


1 


(55) 


and  continues  with  the  equation 


da  Amlf(amm  a)  ^ 


d<5„ 


da 


(56) 


The  first  derivatives  of  the  cofactors  are  generated  by  a  recurrence  which  starts  with 
the  equation 


ddj 

da 


=  0 


(57) 
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and  continues  with  the  equation 

_ _  d&k-l  ,cr,\ 

da  ~ <Lmk  da  °*-‘  *  EfcT  (58) 

until  k  =  m  -  1.  The  second  derivatives  of  the  principal  minors  are  generated  by  a 
recurrence  which  starts  with  the  equation 


and  continues  with  the  equation 

d 0  d^. 


daz 


The  second  derivatives  of  the  cofactors  are  generated  by  a  recurrence  which  starts 
with  the  equation 


and  continues  with  the  equation 


d2At-i 

daz 


until  k  =  m-  1.  The  nth  order  principal  minor  An  is  the  characteristic  determinant 
|A  -  al|. 

The  determinant  and  its  derivatives  are  utilized  by  a  general  routine  for  the 
determination  of  the  real  and  complex  roots  of  the  determinantal  polynomial  p(a). 

The  product  of  the  characteristic  matrix  and  a  characteristic  vector  is  zero  for 
any  characteristic  root.  Sequential  evaluation  of  the  elements  of  the  characteristic 
vector  v^  shows  that  the  ith  element  is  given  by  the  equation 


K  =  -  -  ~  j  E  t 


The  recurrence  for  u*  starts  with  -  1  and  defines  uniquely  the  sequence  of  elements 
for  vM  as  long  as  the  pivotal  element  is  nonzero.  If  the  pivotal  element  vanishes 

then  all  elements  with  index  less  than  i  are  set  equal  to  zero  unless  the  summation 
also  is  zero.  Sequential  evaluation  of  the  elements  of  the  characteristic  vector  v*' 
shows  that  the  jth  element  is  given  by  the  equation 


V?=  -  — —  E  +  l  | 

aJJ+ 1  +  i  ) 


The  recurrence  for  v "  starts  with  =  1  and  defines  uniquely  the  sequence  of  elements 
for  vu  as  long  as  the  pivotal  element  aiJ+1  is  nonzero.  If  the  pivotal  element  vanishes 
then  all  elements  with  index  more  than  j  are  set  equal  to  zero  unless  the  summation 
also  is  zero. 

Two  vectors  are  computed  by  back  substitution.  The  first  element  of  one  vector  is 
unity  and  the  first  element  of  the  other  vector  is  zero.  Each  element  of  the  vectors 
is  the  quotient  of  a  summation  and  a  pivot.  Whenever  the  pivot  is  zero,  the  first  vector 
is  replaced  by  that  linear  combination  of  the  two  vectors  which  reduces  the  summation 
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to  zero,  and  the  second  vector  is  set  equal  to  zero.  The  next  element  then  is  set  equal 
to  zero  in  the  first  vector  and  equal  to  unity  in  the  second  vector.  Finally  the  first 
vector  is  replaced  by  that  linear  combination  of  the  two  vectors  which  satisfies  the 
last  row  of  the  matrix. 

Trianguloidization 

In  the  trianguloidization  of  a  matrix  by  a  similarity  transformation  the  matrix  is 
reduced  one  column  at  a  time.  For  each  successive  diagonal  element,  the  next  element 
above  in  the  same  column  is  used  as  a  pivot  to  eliminate  all  remaining  nonzero 
elements  above.  It  is  required  that  once  elements  are  reduced  to  zero  in  any  cycle 
of  transformation  they  shall  remain  zero  in  all  subsequent  cycles. 

Let  the  partially  transformed  matrix  be  expressed  by  the  equation 


bil  • 

6«  ■ 

•  6U  ■ 

0 

bil  • 

■  b<i  • 

'  64.  ■ 

■  0 

6,1  ' 

O* 

• 

'  6„  • 

0 

bfu  ' 

•  6n,  • 

'  6n.  ■ 

^nn 

where  n  -  s  columns  have  been  trianguloidized.  The  elements  in  the  ith  row  are  replaced 
in  accordance  with  the  transformation 


6<,  b<i  ~  b*-u  O'  >  *  +  1)  (66) 

After  elimination  by  subtraction  of  a  fraction  of  the  (s  -  l)th  row  from  the  ith  row, 
the  same  fraction  of  the  ith  column  is  added  to  the  (s  -  l)th  column.  The  elements 
in  the  (s  -  l)th  column  are  replaced  in  accordance  with  the  transformation 


6*.,-,  -  bk,8—  1 


6»-i,» 


(67) 


The  accuracy  of  the  elements  is  retained  to  within  rounding  error  if  the  elements 
satisfy  the  criteria 


and 


q  s  |6t»lib»-i,>l 
|6yllb«-i.«l 


§  1 


(68) 


n  £  IfrjJIfytil 

6,-1, ,1  bk,»-\ 


§  1 


(69) 


The  pivot  6,_,  ,  is  superior  to  the  pivot  6t,  if  the  pivots  satisfy  the  criterion 

\bi/\  _  16,-i.jl  l6fc,»-il  _  foul 

0  i  V  ib«»l  .  y  !bi.l  lb.-t..l 

|6<jl  +  |6,_i^l  |t>fc,,-il  +  |6*tl 

|6<,l  |6^| 


(70) 


The  similarity  transformation  requires  that  for  any  scaling  of  the  ith  row  there  must 
be  inverse  scaling  of  the  ith  column.  Within  this  limitation  on  scaling  the  pivot  selection 
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is  independent  of  scaling.  \ 

Programming 

SUBROUTINE  RVMTX  (AA,  NA,  RA,  VA.  SA) 

************* *************************** * #*** ♦ ** * #♦ * ************************************** * * * 

FORTRAN  SUBROUTINE  FOR  ROOTS  AND  VECTORS  OF  A  MATRIX 

******************** ^♦♦•••♦•♦••♦••♦•**********************'***********************************  j 

The  address  of  the  matrix  A  is  the  n  x  n  array  AA.  The  order  n  of  the  matrix  is  given  j 

in  the  argument  NA  The  matrix  is  trianguloidized  by  permutations  and  eliminations.  j 

The  real  and  the  imaginary  parts  of  the  characteristic  roots  are  stored  in  alternate  ■ 

addresses  of  the  2n  array  RA  The  covariant  vectors  of  the  matrix  are  stored  columnwise  j 

with  the  real  and  the  imagine  y  parts  of  each  vector  in  alternate  columns  of  the  j 

n  x  2n  array  VA.  The  similarity  ransformation  for  the  trianguloidization  is  stored  in  j 

the  n  x  n  array  SA.  References  < re  made  to  Subroutine  MXRT  for  the  determination  of  i 

roots. 

SUBROUTINE  CHDMT  (MO.  AA,  DA,  1  0,  SD.  AM,  NM,  PM,  DM,  SM)  j 

***«************•******•••*•*«*#•••**•«»»***•*•*••*******************************************  1 

FORTRAN  SUBROUTINE  FOR  CHARACTERISTICS  OF  TRIANGULOID  MATRIX 
*********************************************************************************************  1 


The  mode  of  operation  is  given  in  the  argument  MO.  The  real  and  the  imaginary  parts  j 

of  the  argument  a  are  given  in  the  2-array  AA.  If  MO  is  1,  2,  3  the  computations  are  ; 

carried  as  far  as  the  determinant,  its  first  derivative,  or  its  second  derivative.  The  i 

determinant,  its  first  derivative,  and  its  second  derivative  are  stored  in  the  2-arrays  j 

DA,  DD,  SD.  The  address  of  the  matrix  A  is  the  n  x  n  array  AM.  The  order  of  the  matrix  j 

is  given  in  the  argument  NM.  The  principal  minor  determinants  and  their  first  and  j 

second  derivatives  are  stored  in  the  2n-arrays  PM,  DM.SM.  | 

SUBROUTINE  MXRT  (CD,  CE,  AZ,  NA,  AA,  AM,  PM,  DM,  SM)  ( 

************************* ******************** ************************************************ 

FORTRAN  SUBROUTINE  FOR  CHARACTERISTIC  ROOTS  OF  MATRIX  j 

********************************************************************************************* 


The  step  length  6  for  hunting  is  given  in  argument  CD,  and  the  tolerance  e  for  homing 
is  given  in  argument  CE.  The  initial  position  z  is  given  in  the  2-array  AZ.  The  number 
of  roots  n  is  given  in  argument  NA.  The  address  of  the  matrix  A  is  the  nx  n  array  AM. 
The  addresses  of  principal  minors  and  their  first  and  second  derivatives  are  the 
2n-arrays  PM,  DM,  SM.  The  region  of  a  root  is  located  by  a  hunting  procedure  with 
steps  of  fixed  length,  then  the  root  is  determined  by  a  homing  procedure  with 
Newton-Raphson  iteration.  References  are  made  to  Subroutine  CHDMT  to  obtain  values 
of  the  characteristic  determinant  and  its  first  and  second  derivatives.  The  real  and 
imaginary  parts  of  the  roots  are  stored  in  alternate  addresses  of  the  2n-array  AA. 

SUBROUTINE  OMTXV  (NA,  VL,  VR) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  ORTHONORMAL  VECTORS  OF  MATRIX 
********************************************************************************************* 

The  order  n  of  the  matrix  is  given  in  argument  NA.  The  left  and  the  right  vectors  are 
given  in  the  n  x  2n  arrays  VL  and  VR.  All  vectors  after  a  pivotal  vector  in  one  set  are 
orthogonalized  progressively  with  respect  to  the  pivotal  vector  in  the  other  set,  until 
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each  vector  of  one  set  is  orthogonal  to  the  other  vectors  of  the  other  set.  Finally, 
each  vector  of  one  set  is  normalized  with  respect  to  the  same  vector  of  the  other  set 
in  such  a  way  that  the  largest  elements  in  either  vector  have  the  same  magnitude 
and  the  same  or  opposite  direction. 


DISCUSSION 

The  fact  that  FORTRAN  requires  matrices  to  be  stored  in  transposed  form  is  connected 
with  the  choice  of  indices  for  two-index  arrays.  In  FORTRAN  the  first  index  is  the  most 
rapidly  variable,  whereas  it  would  be  more  natural  for  the  last  index  to  be  the  most 
rapidly  variable.  The  designers  of  FORTRAN  made  the  wrong  choice. 

The  conjugate-gradient  method  was  presented  originally  by  Hestenes  and  Stiefel8 
for  symmetric  positive  definite  matrices.  A  generalization  for  asymmetric  matrices  is 
derived  in  Appendix  A.  Unfortunately,  the  conjugate-gradient  method  is  not  suitable 
for  a  general  purpose  subroutine  because  rounding  errors  cause  rapid  deterioration 
of  the  recurrence  relations. 

Rotations  have  been  used  by  Givens8  and  reflections  have  been  used  by  Householder7 
to  reduce  a  general  matrix  to  triangular  form.  That  Gauss  elimination  is  more  efficient 
for  triangularization  has  been  reported  by  Wilkinson. 

The  method  of  minimized  iterations  was  presented  originally  by  Lanczos8  for 
asymmetric  matrices.  The  algorithm  is  rederived  in  Appendix  B.  Unfortunately,  the 
method  of  minimized  iterations  is  not  suitable  for  a  general  purpose  subroutine 
because  rounding  errors  in  nearly  zero  divisors  upset  the  recurrence  relations. 

Rotations  have  been  used  by  Givens12  and  reflections  have  been  used  by  Householder12 
to  reduce  a  symmetric  matrix  to  tridiagonal  form.  A  generalization  of  the  Householder 
method  to  nonsymmetric  matrices  is  derived  in  Appendix  B.  The  method  fails  if  the 
inner  product  of  a  row  and  a  column  becomes  zero. 

The  figure  of  merit  of  a  matrix  subroutine  is  the  number  of  matrices  in  a  set  of 
test  matrices  which  the  subroutine  is  able  to  process  successfully.  There  was  a  set  of 
test  matrices  on  NORC,  and  there  is  a  set  of  test  matrices  in  a  book  by  Gregory  and 
Karney14.  Matrices  from  both  sets  have  been  collected  into  a  single  set  by  N.  M  Wolcott. 
The  matrices  which  are  defined  for  arbitrary  order  are  cataloged  in  Appendix  C.  To 
these  have  been  added  special  matrices  to  make  a  set  of  fifty  test  matrices  of  the 
sixth  order. 

Extensively  used  in  this  laboratory  is  a  subroutine  MINVR,  which  uses  Gauss-Jordan 
elimination  with  complete  pivot  selection.  The  matrix  must  be  compact  but  is  inverted 
in  place.  The  CDC  6600  subroutine  MATRIX  inverts  a  matrix  in  place  by  Gauss-Jordan 
elimination  with  a  choice  of  no  pivoting,  partial  pivoting,  or  complete  pivoting.  The 
matrix  must  be  compact  but  may  be  a  submatrix  of  a  larger  matrix.  The  IBM  360 
subroutine  MINV  inverts  a  matrix  in  place  by  Gauss-Jordan  elimination  with  complete 
pivoting.  The  matrix  must  be  compact.  The  IMSL  subroutine  UNV1F  uses  the  Crout 
algorithm  to  solve  a  system  of  simultaneous  equations  for  which  the  known  elements 
are  the  columns  of  the  identity  matrix.  The  space  for  two  matrices  is  required,  one 
to  hold  the  original  matrix  and  one  to  hold  the  identity  matrix.  The  matrices  must 
be  compact,  but  may  be  submatrices  of  larger  matrices.  A  subroutine  CROUT  has  been 
prepared  by  A.  H.  Morris,  Jr.  It  uses  the  Crout  algorithm  with  partial  pivoting  to  invert 
a  matrix  in  place.  The  matrix  must  be  compact,  but  may  be  a  submatrix  of  a  larger 
matrix.  The  accumulation  of  inner  products  is  in  double  precision. 

It  is  known  that  for  matrix  inversion  the  Crout  algorithm  requires  the  same  number 
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of  operations  as  the  Gauss-Jordan  algorithm.  The  accumulation  of  each  sum  of  inner 
products  occurs  in  a  single  operation  in  the  Crout  method  and  more  accuracy  is 
possible  in  a  computer  with  a  double-precision  accumulator.  When  the  two  methods 
are  programmed  in  FORTRAN,  their  relative  efficiencies  depend  upon  the  extent  to 
which  the  compiler  takes  advantage  of  parallel  processing  in  the  central  processing 
unit. 

In  a  series  of  test  runs  the  six  subroutines  were  applied  to  the  inversion  of  the 
fifty  test  matrices.  The  subroutine  MATRIX  ran  the  fastest  because  it  is  written  in 
COMPASS.  The  subroutines  CROUT  and  MTRX  were  nearly  equally  fast.  The  accuracy  of 
inversion  was  determined  from  the  product  of  each  matrix  with  its  inverse.  The  extent 
to  which  the  product  deviates  from  the  identity  matrix  is  an  indication  of  the  accuracy 
of  inversion.  All  of  the  subroutines  gave  comparable  accuracy.  CROUT  was  a  fraction 
of  a  digit  better  than  MTRX  for  most  matrices,  but  MTRX  was  better  than  CROUT  for  a 
few  matrices. 

The  CDC  subroutine  MATRIX  uses  the  Householder  method,  and  the  IBM  subroutine 
EIGEN  uses  the  Jacobi  method  to  find  the  roots  and  vectors  of  a  real  symmetric  matrix. 
However,  symmetric  matrices  are  not  within  the  scope  of  the  present  investigation. 
The  IMSL  subroutine  EIGRF  and  the  E1SPACK  subroutines  use  the  QR  algorithm.  A  driver 
EIGV  for  the  EISPACK  subroutines  has  been  prepared  by  A.  H.  Morris,  Jr. 

In  a  series  of  test  runs  the  three  sets  of  subroutines  were  applied  to  the  computation 
of  the  roots  and  vectors  of  the  fifty  test  matrices.  The  subroutines  EIGRF  and  EIGV  ran 
more  than  three  times  as  fast  as  RVMTX.  The  accuracy  of  computation  was  determined 
from  the  difference  between  the  products  of  the  matrix  and  its  vectors  and  the 
products  of  the  vectors  and  their  roots.  With  normalization  of  the  vectors  to  unit 
magnitude,  the  accuracy  of  RVMTX  was  one  digit  better  than  the  accuracy  of  EIGRF  and 
EIGV  for  most  of  the  matrices.  Both  routines  gave  repeated  vectors  for  the  defective 
matrices.  The  roots  and  vectors  are  sorted  in  ascending  order  by  RVMTX,  but  the  roots 
and  vectors  are  not  arranged  in  any  order  by  EIGV. 


CONCLUSION 

Pivot  selection  with  scale  invariance  gives  as  much  accuracy  as  full  pivot  selection 
in  single  precision  subroutines.  The  determination  of  roots  and  vectors  by 
Newton-Raphson  iteration  gives  one  digit  better  accuracy  than  the  determination  by  ■■ 

blockwise  triangularization.  ! 

{ 
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APPENDIX  A 


MATRIX  INVERSION 


< 


CONJUGATE  GRADIENT 


Let  a  circle  be  inscribed  in  a  square  and  let  the  square  and  the  circle  be  deformed 
by  a  linear  transformation  into  a  parallelogram  and  an  ellipse.  The  points  of  tangency 
between  parallelogram  and  ellipse  are  opposite  ends  of  conjugate  diameters.  One 
conjugate  diameter  is  orthogonal  to  the  normal  to  the  ellipse  at  the  end  of  the  other 
conjugate  diameter.  If  the  end  of  the  diameter  and  the  normal  to  the  ellipse  are  given, 
then  the  center  of  the  ellipse  is  in  the  direction  of  the  diameter.  The  ellipse  is  the 
contour  of  constant  value  for  a  quadratic  function  of  the  coordinates.  The  quadratic 
function  is  zero  at  the  center  of  the  ellipse. 

Let  x  be  an  unknown  vector  and  let  y  be  a  known  vector.  Let  r  be  the  residual 
which  is  defined  by  the  equation 


r  =  y  -  A  x  (1) 

The  vector  x  is  a  solution  of  the  equation 

A  x  =  y  (2) 

at  the  point  where  the  quadratic  function  r  ■  r  is  zero.  The  normal  to  the  contour  of 
constant  r  ■  r  is  in  the  direction  of  the  gradient 

V(r  •  r)  =  -  2A'  •  r  (3) 

which  defines  the  direction  conjugate  to  the  vector  r. 

In  the  method  of  Hestenes  and  Stiefel,  the  solution  is  approached  by  a  sequence 
of  vectors  xt,  xn  which  minimize  the  function  r  r.  Associated  with  the  vectors  is  a 
sequence  of  residuals  r,,---,rn  and  a  sequence  of  gradients  A’ ■  r,.  A’ ■  rn.  Let  the 
vector  pt  be  a  linear  combination  of  the  first  k  of  the  vectors  A'  •  •••.  A’  ■  rn.  Let  the 
(k  +  l)th  vector  x**,  be  expressed  in  terms  of  the  fcth  vector  xk  by  the  equation 


I 


=  x*  +  “*Pt  (4) 

where  ak  is  a  scalar.  The  residuals  are  related  in  accordance  with  the  equation 

r*M  =  rk  -  akA  p*  (5) 


The  square  of  the  residual  is  given  by  the  equation 

r»+i  rt+l  =  rk  rk  -  2akpk  A'  rk  +  akzpk  A'  A  pk  (6) 

The  square  of  the  residual  is  a  minimum  where  the  derivative  with  respect  to  ak  is 
zero,  or  where  ak  is  given  by  the  equation 


Pk  A'-rt 
p*  A'  •  A  pt 


(7) 


Substitution  of  this  value  of  ak  in  the  expression  for  rt+]  and  multiplication  by  pt  •  A' 
shows  that 


P»  A’  rtM  =  0 


(8) 


Let  the  vector  p*  be  expressed  by  the  equation 

*-i 

P*  =  A'  •  r*  +  £  fimpm  (9) 

m-t 

where  (3m  is  a  scalar.  Let  the  vectors  be  orthogonal  as  expressed  by  the  equation 

Pn  A'  A  p,  =  0  (m  *  k)  (10) 


Then  pm  is  given  by  the  equation 


Pm  A'  A  A'  rt 
Pm  ■  A*  •  A  ■  pm 


(11) 


and  the  vectors  are  constructed  by  a  multi-term  recurrence. 

For  any  m<k  the  residual  rt+1  can  be  expanded  in  the  polynomial 

r*+l  =  rm+l  —  am  +  lA  ■  Pm+1  —  ~  OtA  •  p*  (12) 

Multiplication  by  pm  •  A'  leads  to  the  equation 

Pm  A’  rt+1  =  0  (13) 

Thus  the  square  of  the  residual  is  a  minimum  with  respect  to  any  fix  given  by  the 
equation 

* 

<5X  =  E  Pm«5«m  (14) 

m=l 

The  square  of  the  residual  is  reduced  to  zero  after  n  steps. 

For  any  m  <  k  the  vectors  are  related  by  the  equation 

r*  A  pm  =  0  (m  <  k)  (15) 

or  for  m  =  Jfc  by  the  equation 

r*  A  p*  =  r*  A  A'  r*  (m=fc)(16) 

and  for  m  #  k  by  the  equation 

r»  •  A  •  A'  rm  =  rm  •  A  A’  •  r*  =  0  (m*k)  (17) 


Substitution  for  amA  ■  pm  its  expression  in  terms  of  residuals  leads  to  the  equation 


(fmu  -  I'm)  •  A  •  A'  •  rk 

(rm+l  —  ^m)  '  A  '  Pm 


(18) 


from  which  it  is  apparent  that  fim  =  0  for  all  mcfc-l.  Thus  the  vectors  pt  are 
synthesized  by  a  two- term  recurrence  equation 

P*  =  A'rt  +  Pi-iPk-i  (19) 


for  which  is  given  by  the  equation 

_  r*  A  A'  rk 
P*-‘  ~  r*_,  A  A'  r*_, 


(20) 


The  effect  of  rounding  error  on  the  mutual  orthogonality  of  the  vectors  tends  to  grow 
rapidly  with  the  two-term  recurrence,  but  can  be  checked  with  the  multi-term 
recurrence. 

Inasmuch  as  the  vectors  A  p,  are  mutually  orthogonal,  they  can  be  used  to  express 
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,  t 

- 


the  identity  matrix  by  the  equation 


I  =  f  P»  A  A  P» 
*.»  P»  A'  A  •  p* 

whence  the  inverse  matrix  is  given  by  the  equation 

A-i=  y  P*A  P* 
ti  P*  A'  A  •  p4 


(21) 

(22) 


The  accuracy  of  the  inverse  matrix  is  limited  by  the  extent  to  which  the  initial  residual 
vector  has  substantial  components  in  the  directions  of  all  of  the  characteristic  vectors 
of  the  matrix. 
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APPENDIX  B 


CHARACTERISTIC 

ROOTS  AND  VECTORS 


COMPANION  MATRIX 


A  companion  matrix  has  ones  on  the  diagonal  next  above  the  main  diagonal,  nonzero 
elements  along  the  bottom  row,  and  zero  elements  elsewhere.  The  characteristic 
equation  of  the  matrix  A  is 

(-1)"|A  —  a  I|  =p(a)  =  £  cmam  (1) 

m=0 

where  the  nth  order  coefficient  of  the  characteristic  polynomial  p(a)  is  unity.  The 
characteristic  equation  of  a  companion  matrix  is 


-a 

1 

0 

0 

0 

0 

0 

-a 

1 

0 

0 

0 

0 

0 

-a 

0 

0 

0 

0 

0 

0 

-a 

1 

0 

0 

0 

0 

0 

-a 

1 

~c0 

-Cl 

-e2 

-c„-3 

~Cn-2 

-c„-,~a 

That  the  characteristic  determinant  is  equal  to  (-l)”p(a)  may  be  verified  if  each 
column  on  the  right  is  multiplied  by  a  and  is  added  to  the  next  column  on  the  left 
until  diagonal  elements  above  the  bottom  row  are  eliminated,  and  columns  are 
interchanged  to  bring  the  first  column  into  the  position  of  the  last  column.  Then  the 
matrix  of  the  determinant  is  triangular  and  the  determinant  is  equal  to  the  product 
of  its  diagonal  elements,  which  are  all  ones  except  the  last  element. 

The  product  of  the  characteristic  matrix  and  a  characteristic  vector  is  zero  for 
any  characteristic  root.  Sequential  evaluation  of  the  elements  of  the  characteristic 
vector  vM  shows  that  the  ith  element  is  given  by  the  equation 


1 


a 


t-i 


(3) 


a 


n- 1 


Sequential  evaluation  of  the  elements  of  the  characteristic  vector  vv  shows  that  the  j th 
element  is  given  by  the  equation 

v*  =  IIPn-i(a„).  ?„-,(<*„).  •••.  HI  (4) 

where  p*(a1/)  is  just  the  fcth  order  polynomial  which  is  constructed  during  the  nested 
method  of  polynomial  evaluation. 
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If  all  characteristic  roots  are  distinct,  then  the  vectors  are  proportional  to  the 
columns  of  a  modal  matrix  V,  and  the  vectors  vv  are  proportional  to  the  rows  of  the 
inverse  matrix  V-1. 


REDUCTION  TO  COMPANION 

In  the  Danilevsky  reduction  a  matrix  is  transformed  into  its  companion  matrix  by 
a  sequence  of  eliminations.  The  first  element  above  each  diagonal  element  is  used  as 
a  pivot  to  eliminate  all  other  elements  in  the  same  row.  Let  the  partially  transformed 
matrix  be  expressed  by  the  equation 


B  = 


0 

0 

0 

0 

0 

1 

0 

0 

b.l 

■ 

b«,«+i 

bwn 

ba*l.l 

b«+i,» 

bt+»,«+i 

^•  +  l.n 

bnl 

bn. 

bn.»  +  l 

t>nn 

(5) 


where  s  -  1  rows  of  the  matrix  have  been  reduced  to  upper  diagonal  form.  The  elements 
in  the  jth  column  are  replaced  in  accordance  with  the  transformation 


■*  bu 


(6) 


After  elimination  by  subtraction  of  a  fraction  of  the  (s  +  l)th  column  from  the  jth 
column,  the  same  fraction  of  the  j>th  row  is  added  to  the  (s  +  l)th  row.  The  elements 
in  the  (s  +  l)th  row  are  replaced  in  accordance  with  the  transformation 


b»*l.t  ■*  b.+  i .*  + 


(7) 


After  elimination  the  (s  +  l)th  column  is  divided  by  6„+,  and  the  (s  +  l)th  row  is 
multiplied  by  6s,,+l- 

The  method  of  reduction  fails  for  a  diagonal  matrix. 

Even  after  the  characteristic  polynomial  has  been  derived  by  reduction  of  the 
matrix,  a  polynomial  of  high  degree  cannot  be  evaluated  for  large  arguments,  and  the 
largest  roots  are  indeterminate. 


TRIDIAGONALIZATION 

In  the  tridiagonalization  of  a  matrix  by  a  similarity  transformation  the  matrix  is 
tridiagonalized  progressively  one  row  and  one  column  at  a  time.  At  each  successive 
diagonal  element,  the  next  right  element  in  the  same  row,  and  the  next  lower  element 
in  the  same  column  are  used  as  pivots  to  eliminate  the  remaining  nonzero  elements 
in  the  row  and  the  column.  It  is  required  that  once  elements  are  reduced  to  zero  in 
any  cycle  of  transformation  they  shall  remain  zero  in  all  subsequent  cycles. 


Let  the  partially  transformed  matrix  be  expressed  by  the  equation 


bn 

0 

0  0 

0 

0 

b.. 

b.,.+i  b.j 

b»n 

0 

b»+i.» 

bi+i.«+i 

0 

b4. 

1  0 

bn. 

bnn 

(8) 


where  s  -  1  rows  and  columns  have  been  tridiagonalized 

The  Householder  tndiagonalization  is  based  on  the  application  of  a  similarity 
transformation  in  which  the  prefactor  and  the  postfactor  both  are  of  the  same  form, 


I  -  2 


uv 
U  V 


The  prefactor  and  the  postfactor  satisfy  the  identity 


;  uv  \  /  uv  \ 

1-2 -  !~2 -  =1 

V  uv/\  uv/ 


(9) 


(10) 


The  application  of  this  similarity  transformation  to  matrix  B  is  expressed  by  the 
equation 


(l  -  2  — — •  B  •  /i  -  2  — -V  ) 

\  u  V  /  \  u  V  ' 


2  **(v  •  B)  _  2  (B  u)v 


U  V 


U  V 


(U) 

(u  v)* 


If  all  elements  of  index  less  than  s  -  1  are  zero  in  both  u  and  v,  then  uv  contains  no 
elements  in  the  sth  column  or  sth  row,  but  u(v  B)  contains  elements  in  the  sth 
column,  while  (B  u)v  contains  elements  in  the  sth  row  For  elements  below  and  to 
the  right  of  the  pivotal  elements  to  become  and  remain  zero,  it  is  necessary  for  the 
elements  of  u  and  v  beyond  the  pivotal  element  to  be  proportional  to  the  column  or 
the  row  which  is  to  be  eliminated. 

Let  the  vector  u  be  the  column  vector 


I  0 

r 

0 

u  -  ,  (12) 

b.*2,> 


6 


m 


and  let  the  vector  v  be  the  row  veetor 

v  9.  ,  0,  *  ^*n 


(13) 


1 


* 


3 


(14) 


The  scalar  product  of  u  and  v  is  given  by  the  equation 

TV 

«  v  =  <j.,.+1q,+ll.  +  E  b,kbkt 
*>=•+2 

The  ith  element  of  the  sth  column  becomes 

n 

<7«.»+ib,+i.,  +  E  btkbk, 
bu-2 - -  6*.  (i*  s  +  1)  (15) 

9*.»+t9«+i,*  +  E  b,kbkt 
*=•+2 

Thus  the  elements  in  the  sth  column  are  eliminated  if  the  parameters  9,.,+1  and  q.+ 
satisfy  the  equation 

TV 

9*,#+i9*+t,»  —  1  b>  +  1  #  —  ^  —  0  (16) 

4  =  #  4-2 

The  jth  element  of  the  sth  row  becomes 


6s>-  2 - *=£* -  0'**+  l)  07) 

<7c,a+l9»  +  l.a  +  S  ^afc^ka 

fc=a  +  2 

Thus  the  elements  in  the  sth  row  are  eliminated  if  the  parameters  q.,,M  and  q,+1, 
satisfy  the  equation 

TV 

9*.**i<?**i,»  —  26,,,  +  i9,+  i,s  -  E  b,kbk,  --  0  (18) 

t-a+2 

Comparison  of  the  equations  for  the  parameters  shows  that  they  are  reduced  to  a 
single  equation  by  the  substitutions 

9**1, a  =  Vb,+  U,  9*. *41  “  vbs,a+l  (19) 

where  the  parameter  tj  is  a  solution  of  the  equation 


V  2  -  2V 


E  b,kbk, 

k=t* 2 _ 

^*,■*1^*41,* 


=  0 


Thus  t)  is  given  by  the  equation 


q  =  1  ± 


E  bwkbk. 


b»,»*ib*4i,. 


1 

2 


(20) 


(21) 


and  u  ■  v  is  given  by  the  equation 

i 

u  v  =  ±  8(6,i.+16ttl.,)l(  E  2  E  b,kbk,  (22) 

't-*4l  /  fl<l 

If  the  product  of  the  pivotal  elements  happens  to  be  zero,  then  rows  and  columns 


t  t 
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can  be  permuted  to  bring  into  pivotal  position  a  pair  of  elements  with  a  finite  product. 
The  value  of  T)  is  real  only  if  the  product  6,,,+i6,+i,,  has  the  same  sign  as  the  sum 

E  b,A.  (23) 

t-«+i 

There  always  must  be  at  least  one  pair  of  elements  for  which  the  product  of  the  pair 
has  the  same  sign  as  the  sum  of  products.  If  there  is  more  than  one  the  pair  with 
the  product  of  largest  magnitude  is  selected  for  pivotal  elements.  The  sign  of  the 
radical  is  selected  to  give  u  v  the  larger  magnitude. 

The  method  fails  if  the  sum  of  the  products  is  zero  as  expressed  by  the  equation 

E  b,kbk.  =  0  (24) 

Then  no  permutations  or  eliminations  can  prevent  u  v  from  having  a  value  of  zero. 
The  method  fails  for  a  companion  matrix,  which  may  have  distinct  roots  and  vectors. 
The  companion  matrix  contains  information  only  about  roots  and  not  about  the  vectors 
of  any  matrix  from  which  the  companion  matrix  may  have  been  derived. 


MINIMIZED  ITERATION 


In  the  method  of  minimized  iteration,  the  roots  and  vectors  are  determined  from 
a  sequence  of  biorthogonal  vectors  in  which  each  successive  vector  is  derived  from 
previous  vectors  through  multiplication  by  the  matrix  A. 

Let  ja,  ■■  ,  jn_,  and  j°,  — ,  jB_1  be  biorthogonal  vectors  such  that 


— i_(. 

r  a  l_,  \ 


A  jM-i 


M-l 

E 

v-0 


}v) 


V 


Iv  J 


(25) 


and 


Multiplication  by  jx  is  given  by  the  equation 


J  j“‘‘  A  j„\J  j7o  j / 

given  by  the  equation 

,x  j  _  __J__  /;>  .  A  j  VjX  A 

J  r-A  v.  V  x  ^  r  / 


(26) 


(27) 


and  multiplication  by  jx  is  given  by  the  equation 


_  1  /,-■  ,  .  yr  ■*  jj-  j>\ 

'  h  r  *  h  i.  i.  r  I 


(28) 


Thus  the  vectors  are  biorthogonal  in  accordance  with  the  equations 

f  j„  =  0  J"  jx  =  0  (29) 

for  any  A  <  /U  if  jx  •  j„  =  0  or  jv  •  jx  =  0  for  v  *  A.  Since  the  vectors  are  biorthogonal  for 
p.  =  1,  they  are  biorthogonal  for  any  n  by  induction.  The  vectors  are  normalized  in 
accordance  with  the  equation 

r  =  i  oo) 
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Multiplication  by  jx  ‘  A  gives  the  equation 

jX~‘  A-  }M  =  —  (jx"‘  •  A'  A-  j„_,  -  E  jx'1 '  A  jj*'  A  •  jM_,) 

and  multiplication  by  A  ■  gives  the  equation 

jX  ■  A  j„_,  =  .x-f1-  -'  (j’1’1  A  A  jM_,  -  E  jX'1  A  jj*  ■  A  ■  j„_,) 

J  A  Jx  \  p-=0  / 

Elimination  of  jx_1  •  A  -  A  •  leads  to  the  equation 

£  (^AjJff  A-v,)^ 

l/*=X+l 

This  equation  can  be  true  for  all  X  g  fj.  -  1  only  if  the  vectors  satisfy  the  equation 

jx  A  j„  =  0  (X<M  -  1)  (34) 

Multiplication  by  j'1"1  A  gives  the  equation 


(31) 


(32) 


(33) 


jM-i  A  jx  =  x  /  ;  (j"  1  '  A  A  jx_!  -  E  jM“‘  A  jj1'  ■  A  jx-,) 

J  ‘  A  ‘  Jx-i  \  i/=o  / 


and  multiplication  by  A  jX- 1  gives  the  equation 


r  a  i-,  -  (, 


A  A  jx_,  -  E  j"'1  A  jj1'  A  jx_, 


(35) 


(36) 


Elimination  of  jM  1  A  A  jx.,  leads  to  the  equation 

E  (j“"‘  4  A  j„)(j1'  •  A  jx_,)  =  0  (37) 

l>  =  \+t 

This  equation  can  be  true  for  all  \  1  only  if  the  vectors  satisfy  the  equation 

J*  A  jx  =  0  (\</i-l)(38) 

Thus  the  vectors  satisfy  the  three-term  recurrence  equations 


Jm 


j<T  /  j  —  (a-V.  -  j„-ijM"1  A  j„_,  -  j^-zj'*'2-  A  • jM_,) 


and 


j"  =  pi-",-  (j""1  A  -  j""1  A  j^J""1  -  j"'1  A  jM_,j'*-*) 
Let  (?M_!  be  defined  by  the  equation 

1  =  j"'*  ■  A  jM_, 

and  let  be  defined  by  the  equation 

7„-i  =  jM  A  •  j„.,  =  jM_l  A  jM 
Then  the  three-term  recurrence  equations  become 


j<*  =  7m-i 


Pn-iin-i  7u 


(39) 

(40) 

(41) 

(42) 

(43) 
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and 


During  each  cycle  of  iteration  the  value  of  is  computed  from  the  square  root  of 
the  scalar  product  of  the  unnormalized  vectors. 

Insofar  as  the  nth  vectors  j„  and  jn  are  orthogonal  to  n  independent  vectors,  they 
must  be  zero  within  rounding  error.  Inasmuch  as  each  vector  is  obtained  from  j0  and  j° 
by  powers  of  the  matrix  A,  the  nth  vectors  must  be  the  products  of  the  vectors  j0  or  j° 
and  the  characteristic  polynomial  p( A)  because  by  the  Cayley-Hamilton  theorem  the 
matrix  A  is  a  solution  of  its  .own  characteristic  equation. 

Let  a  sequence  of  polynomials  p0(a),  ■■■ ,  pn(a)  be  generated  by  a  recurrence  which 
starts  with  the  equation 


Po(«)  =  1 


and  continues  with  the  equation 


P/i(a)  =  - -  i(«)  “  ^-iP^-i(a)  -  TV-aP^-ata)  (46) 

7(i-i  t  ) 

until  pn(a)  =p(a).  The  nth  polynomial  and  its  derivatives  can  be  utilized  by  a  general 
routine  for  the  determination  of  the  real  and  complex  roots  of  the  characteristic 
polynomial. 

Insofar  as  the  matrix  A  is  expressed  by  the  equation 

n 

A  =  £  (47) 

t-=i 

a  polynomial  p^( A)  is  expressed  by  the  equation 

n 

Pi i(A)  =  £  P„(a„) a^a-  (48) 

V=l 

The  vectors  j„  and  are  expressed  by  the  equations 

Jd  =  Pd(A)  jo  j“  =  j°  PM( A)  (49) 

Inasmuch  as  the  vectors  and  jM  are  biorthogonal,  they  can  be  used  in  an  expansion 
of  the  identity  matrix  I  in  accordance  with  the  equations 


I  =  £ 


I  =  £  nv 


Multiplication  of  the  characteristic  vectors  by  the  identity  matrix  and  substitution  for 
vectors  their  expressions  in  terms  of  polynomials  lead  to  the  equations 

n 

a(i  =  (j°  a(i)  £  Pv(an)iv  (5i) 


aM  =  (jo  '  aM)  £  Pw(«d)j1'  (52) 

v=i 

The  accuracy  of  the  roots  and  vectors  is  limited  by  the  extent  to  which  the  initial 


t 


vectors  have  substantial  components  in  the  directions  of  all  of  the  characteristic 
vectors  of  the  matrix. 

If  the  matrix  has  a  multiple  root,  it  is  necessary  to  repeat  the  iteration  with  new 
initial  vectors  until  all  of  the  vectors  have  been  determined  for  the  multiple  root. 


APPENDIX  C 


TEST  MATRICES 


CATALOG  OF  MATRICES 
Test  Matrices  for  n  =  6  in  Card  Decks 

Let  A  =  identity  matrix. 

a*,  =  0  (i  j) 

=  1  (i  =  j) 

A,  =  A  =  I  =  A'1 

Let  A  =  circulative  permutation. 

atj  =  0 
ai}  =  1 

A2  =  A 

Let  A  =  reversive  permutation. 

=  0 
aij  ~  1 

A3  =  A  =  A'  =  A-1 

Let  A  =  asymmetric  permutation. 

aiy  =  0  (j  *  (n  -  Zi  +  l)modulo(n  +  1)) 

a4,  =  1  (J  =  (n  -  Zi  +  l)modulo(n  +  1)) 

A4  =  A 

Let  A  =  positive  definite  circulative  matrix. 


(j  *  n  +  i  -  i) 
(j  =  n  -!  1  -  i) 


(j  *  (i  +  l)modulo  n) 
(j  =  (i  +  l)modulon) 

A-1  =  A’ 


Let  A  =  inverse  circulative  matrix. 


ai>  -  - 

=  + 


1 

2  n* 

2n  -  1 
2n2 


A6  =  (72)A  =  (72)Aj 1 


Let  A  =  tridiagonal  finite  difference  matrix. 


a\j  =  0 

aH  =  -  1 


A7  =  A 


Let  A  =  inverse  tridiagonal  matrix. 


n  +  1  -  j 

a-u  = - —  i 

’  n  +  1 


n  +  1  -  i 

ati  = - j 

*  n  +  1  J 


A8  =  (7)A  =  (7)Af 1 

Let  A  =  symmetric  Frank  matrix. 


a*,  =  min(i,;) 


A,  =  A 


Let  A  =  inverse  Frank  matrix. 

a«;  =  0 

at,  -  1 

S'2 
aw  -  +■  1 

Ajo  --  A  —  Ag 


O'  *  i) 

0  =  i) 


(l#-i|>  1) 
i) 

(l?  -  i|  =  0) 


0  ^  *) 
(i  s ;') 


(j  i  n 
(j  i  n 

(i  j  '  n) 
{i  j  n) 


2 


..  — KM***  . 


Let  A  symmetric  Boothroyd  matrix 

ai}  =  max(i.j) 


A,,  A 


Let  A  inverse  Boothroyd  matrix 


0 

A 

1 

ai)  =-- 

+  1 

(b-il=  1) 

aij  ~ 

-  1 

(i  =j=  1 ) 

a<,  = 

+  2 

(x  =  j,  1  <  x  <  n) 

“i I  = 

n  -  1 

n 

«*>. 

II 

II 

A12  (6)A  -  (6)AU* 

Let  A  symmetric  Givens  matrix. 


A,:,  A 


ai}  ~  2  min(x,;)  -•  1 


Let  A  inverse  Givens  matrix. 


av  = 

0 

(iy-<i  >  i) 

ai>  = 

1 

z 

(b  -  il  =  1) 

aiJ  = 

+  § 

(x  =  ;  =  1 ) 

aij  ^ 

+  i 

(i  =  j,  l  <  i  <  n) 

a4j  = 

+  1 

(i=3  =  n) 

A„  (2)A  (2)A13‘ 


3 


Let  A  =  symmetric  Lehmer  matrix. 


* 

1 

a«  =  J 

VII 

• 

i 

-  • 

0  s  i ) 

A.»  =  (60)A 

Let  A  =  inverse  Lehmer  matrix. 

=  0 

(1;  - 1|  >  l) 

i  +  1 

*  2i  +  1 

O'  =  i  +  1 ) 

j+  1  . 

"  2;  +  1  ; 

0  =  4-1) 

4i3 

w  4i2  -  1 

(t  ~  j.i<  n) 

n2 

i}  2n  -  1 

(t  =  ;  =  n) 

A16  =  (3465)A  =  (207900)A1"51 

Let  A  =  symmetric  Todd  matrix. 

ati  =  |  i-  j\ 

A17  =  A 

Let  A  =  inverse  Todd  matrix. 

(1  <(>  —  i|<n  —  1) 
(1 i  -  tl  =  i) 

-  i|  =  n  -  1) 

(i  =;  =  1) 

(i=j.  1  <  i  <  n) 

(i  =  j  -  n) 

A,«  =  (10)A  =  (10)Ar7' 

4 


a4,  =  0 


atj  =  + 


1 

_e _ 

n  -  1 

n  -  2 
2n^2 


ay  =  -  1 


n  -  2 
2n  -  2 


I 


Let  A  -  symmetric-  Lietzke  matrix 


«*>  =  n  [i  j 


Aio  A 


Let  A  inverse  Lietzke  matrix. 


a4)  =  0 


a«  -  +  ;r- 

x’  2n  +  2 

n  +  2 

a,,  -  +  -  — 

’  2  n  +  2 


a»)  =  +  1 


n  +  2 
a*’  +  2n  +  2 


( 1  <  [7  -  i|  <  n  -  1 ) 


([?  i|  =  1) 


([>  -  t|  =  n  -  1) 


(i  =  j=  1 ) 


(i  =  j.  1  <  i  <  n) 


(i-j  =  «) 


A20  =  (14)A  1  =  (14)A,V 


Let  A  =  symmetric  Pascal  matrix. 


(i  +  j  ~  2)! 
(t-  1)'0 -O' 


A„  =  A*1 


Let  A  -  trigonometric  orthogonal  matrix 


a^(nhYSin(^Ll) 


A,,  A  A  A  1 


Let  A  =  Hilbert  matrix. 


**  "  i  +  j  -  1 


A2t  =  (27720)A 


Let  A  -  inverse  Hilbert  matrix. 


+  l)!(rt  +  j  -  1)! 
(i  +  j-  l)(i  -  1)!*0  -  !)!*(*.  - 


A25  =  A  =  (27720)AaV 


Let  A  =  symmetric  Hadamard  matrix. 

Let  k  be  any  integer  and  let  p  =  n  +  1  where  p  is  a  prime  number.  The  Legendre -Jacobi 
reciprocity  function. 

-  1  if  k  *  mp,  k  *  m2modulop 
Q)='  0  if  A  =  mp 

+  1  if  Ac  =  m2modulop 


where  m  is  any  integer.  Then 


A2«  —  A 


Let  A  =  inverse  Hadamard  matrix. 


1 

/  *  \ 

/  j  'll 

• 

a“  "  n  +  1 

An+  1/ 

U+ 1/ 

\n+  1/J 

A„  =  (7)A  =  (7)A2V 


■ -M:'.  - ****** :  ,-,l,  - 


Let  A  =  random  Forsythe  matrix. 

Oy  =  random  number  m  in  range  -  100  <  m  <  +  100 


^28  —  A 

Let  A  =  Clement  matrix. 


aii  =  0 
dfj  =  i 

aij  =  n  -  j 


(l?-i|*  1) 
0  =  i  +  1) 
0  =  i  -  1) 


Azs  -  A 


Ajo  ~  (15)A2S1 


Let  A  =  linear  circulative  matrix. 

Oy  =  0  -  i  +  1  )modulo  n 

A;ii  =  A 

Let  A  =  inverse  circulative  matrix. 

(i  *  j,  j  *  (i  +  1  )modulo  n) 
0  =  (x  +  l)modulo  n) 

(i  =  ;) 


=  + 


a« 


=  + 


n*(n  +1) 
n2  +  n  +  2 


n2(n  +  1 ) 

(n  -  l)(n  +  2) 
n2(n  +  1 ) 


A32  =  (126)  A  =  (126)A3,‘ 


Let  A  =  permuted  Wilkinson  matrix. 


^  33  ~  A 


0  £  i) 

0'  =  i  +  1) 
0'  >  i  +  1 ) 

A34  =  A'1 


Oy  =  (-ir> 

ai)  ~  1 

Oif  =  0 
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*.  ■* 


Let  A  -  binomial  matrix 


» 


A3s  -  A  A  1 


_  (-l)*-»(t-  1)' 
ax>  (i-i)'(j-l)' 


Let  A  Vandermonde  matrix. 


a 


«  =r 


A»  -  A 


A37  - 


Let  A  =  interpolative  matrix 


A3b  •  A  Aj9  - 

Let  A  -  nondefective  companion  matrix. 

ait  =  0  (i  *  n. 

ait  =  +  1  0  =  * 

ati  =  -  1  (i  =  n' 


A«o  ~  A 


Let  A  =  special  symmetric  Martin  matrix. 
This  matrix  has  three  pairs  of  equal  roots. 
A41  =  A 


Let  A  =  special  symmetric  Voigt  matrix 
This  matrix  has  three  pairs  of  equal  roots 
A«j  =  A 


(120)A„' 


(24)ASV 


;  *  t  +  1) 
i- 1) 


a 


Let  A  =  special  nonsymmetric  Varah  matrix. 

The  roots  are  (-3. -2,3.  1,|.2)  for  which  there  are  six  independent  pairs  of  vectors. 
A*3  =  (6)A 

Let  A  =  defective  companion  matrix. 

at<  =  0  (i  *  n.  j  *  i  +  1) 


atj  =  1 

(-  l)n_,n( 

a<^  {n-m 


0  =  1+1) 


(i  =  n) 


This  matrix  has  six  roots  equal  to  1.  for  which  there  is  only  one  pair  of  vectors. 

A44  =  A  A45  =  A 

Let  A  =  defective  triangular  matrix. 

=  (tij) 

da  =  0  (j>i) 

This  matrix  has  six  roots  equal  to  1,  for  which  there  is  only  one  pair  of  vectors. 

A4e  =  A  A47  =  A 

Let  A  =  special  defective  Voigt  matrix. 

The  roots  are  (-1.  -1.  -1.  -i,  +i,  +1)  for  which  the  three  roots  equal  to  - 1  have  only 
one  pair  of  vectors. 

A4a  =  A 

Let  A  ^  special  defective  Varah  matrix 

The  roots  are  (1.  1,  2  -  i.  2  +  t.  3.  3)  for  which  the  two  roots  equal  to  +1  have  only 
one  pair  of  vectors 

A49  =  A 


Let  A  =  symmetric  singular  matrix. 

The  roots  are  (0,0.  3.4,4,  13)  for  which  there  are  six  independent  pairs  of  vectors. 
Am  =  A 
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